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' Abstract In some fields such as Mathematics Mechanization, automated reasoning and Trustworthy 
, Computing etc., exact results are needed. Symbolic computations are used to obtain the exact results. 
^ ' Symbolic computations are of high complexity. In order to improve the situation, exactly interpolat- 
O . ing methods are often proposed for the exact results and approximate interpolating methods for the 
approximate ones. In this paper, we study how to obtain exact interpolation polynomial with rational 
, coefficients by approximate interpolating methods. 

^ [ Key words Numerical approximate computation, symbolic-numerical computation, continued frac- 
tion, multivariate interpolation, Vandermonde determinant. 

^ : 

^ J i 1 Introduction 
G\ • 

! Some fields such as automated reasoning and trustworthy computing etc., need exact results, 

and symbolic computations are used to obtain the exact results. Symbolic computations are 
principally exact and stable. However, they have the disadvantage of intermediate expression 
swell. Numerical computations have not the problem, however only give approximate results. 
In recent two decades, numerical methods are applied in the field of symbolic computations. In 
; ^ ' 1985, Kaltofen presented an algorithm for performing the absolute irreducible factorization, and 
^ ' suggested to perform his algorithm by floating-point numbers, then the factor obtained is an ap- 
proximate one. After then, numerical methods have been studied to get approximate factors of 
a polynomial^^^^l. In the meantime, numerical methods are applied to get approximate greatest 



Yong FENG 

School of Computer Science and Engineering, University of Electronic Science and Technology of 
China, Chengdu 610054, China; Laboratory for Automated Reasoning and Programming, Chengdu Institute 
of Computer Applications, Chinese Academy of Sciences, Chengdu 610041, C/ima.Email:yongfeng@casit. ac.cn. 
Jingzhong ZHANG 

Laboratory for Automated Reasoning and Programming, Chengdu Institute of Computer Applications, Chinese 
Academy of Sciences, Chengdu dlOOAl, China. 
Xiaolin QIN 

Laboratory for Automated Reasoning and Programming, Chengdu Institute of Computer Applications, Chinese 
Academy of Sciences, Chengdu dlOOAl, Cfcina.Email: qinxl811028@163.com. 
Xun YUAN 

Laboratory for Automated Reasoning and Programming, Chengdu Institute of Computer Applications, Chinese 
Academy of Sciences, Chengdu 6100AI, China. 

*This research is supported by China 973 Project NKBRPC-2004CB318003, the Knowledge Innovation Pro- 
gram of the Chinese Academy of Sciences KJCX2-YW-S02, and the National Natural Science Foundation of 
China(Grant NO. 10771205). 



2 



YONG FENG • JINGZHONG ZHANG • XIAOLIN QIN • XUN YUAN 



common divisors of approximate polynomials'^" to compute functional decompositions!^ ^1, 
to test primalityl^^l and to find zeroes of a polynomial!^'''!, jn 2000, Corless et al. applied 
numerical methods in implicitization of parametric curves, surfaces and hypersurfaces'^^!, and 
the resulting implicit equation is still an approximate one. In [15],Cheze et al. discussed how to 
obtain an exact absolute polynomial factorization from its approximate one, which only involves 
recovering an integral coefficient from its approximation. 

Interpolation methods as an efficient numerical method have been proverbially used to 
compute resultants and determinants, etc. And approximate interpolation methods arc still 
used to get the approximate ones'^^"^*!. In order to obtain exact results, people usually use 
exact interpolation methods to meliorate intermediate expression swell problem arising from 
symbolic computations [i^'^^s] -^u fact, these arc not approximate numerical computations but 
big number computations, which are also exact computations and only improve intermediate 
expression swell problem. Recently, Zhang et al. proposed an algorithm to recover the exact 
rational number from its approximation!^^! , and built a bridge by which exact results can be 
obtained by numerical approximate computations. In this paper, we discuss how the errors 
of support points affect that of coefficients of the interpolating polynomial, thereby present 
an algorithm to use approximate methods to get exact interpolation multivariate polynomial. 
The algorithm can be carried out in parallel computers. Compared with exact interpolation 
method, the efficiency of the our algorithm is higher when the scale of problem is larger. This 
paper provide a way to obtain exact results by numerical computations. 

The remainder of the paper is organized as follows. Section 2 gives a review of modified 
continued fraction method, by which an exact rational number can be obtained from its ap- 
proximation, and gives a review of Kronecker product of two matrix. Section 3 first proposes 
estimation of the error to ensure the exact interpolation polynomial to be obtained, and then 
presents an algorithm to obtain an exact polynomial from its approximation for univariate and 
multivariate polynomial over rational number field, respectively. Section 4 gives some experi- 
mental results. The final section makes conclusions. 

2 Preliminaries 

In general, obtaining exact polynomial by approximate computations consists of two steps. 
First compute approximate polynomial within an error control by approximate numerical meth- 
ods, and then recover the exact coefficients of the approximate polynomial by continued fraction 
method. In this section, we review the main results to be used in this paper. 

A continued fraction representation of a real number x is one of the forms: 

ao H , (1) 

ai H J 

02 H ; 

as H 

where ao is an integer and Oi,a2,a3, • • • are positive integers. One can abbreviate the above 
continued fraction as [ao; ai, a2, • • • ]. Truncating the above continued fraction representation of 
a number x early yields a rational number which is in a certain sense the best possible rational 
approximation. We call [ao;ai,-- - ,a„] the n-th convergent of [ao; ai, 02, • ■ • ]. By continued 
fraction representation, the relation between a rational number and its approximation was 
disclosed as follows 

Theorem 1 Let no/ni be a reduced rational number and r its approximation. Assume that 
no,ni are positive integers and N > max{ni,2}. We have the continued fraction representation 
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r = [bo,bi,--- ,b]\i] and |r — no/rii| < 1/(2A^^). Then, no/ni = [bo,bi,--- ,6^], where either 
M = L or M > L. In the case of M > L, for any positive integer L < t < M , the denominator 
of rational number g = [60, bi, - ■ ■ ,64] is greater than N . 

Based on the theorem above, an algorithm for obtaining the exact number was designed as 



Algorithm 1 Input: a nonnegative floating-point number r and a positive number N ; 
Output: a rational number b. 

Step 1: Set i = Q, = 1, /i_2 ~ 0, fc_i = 0, fc_2 = 1; 'f^c' convert r to a rational 
number and assign it to tern; 

Step 2: Get integral part of tern and assign it to a, assign its remains to b; 

Step 3: Compute hi — a* + ki — a* fci_i + ki^2- If ki > N , then goto 

Step 6; 

Step 4-' Set i :^ i + 1; 

Step 5: If b 0, set tern = ^ and goto Step 2; 
Step 6: Compute hi-i/ki^i and assign it to b; 
Step 7; return b. 

It follows from theorem [T] that once an bound N on the denominator of a positive rational 
number is estimated, we can obtain the exact rational number as follows. Compute its ap- 
proximation with the error less than 1/{2N^), and then recover the rational number from its 
approximation by algorithm [T] 

Now, let us give a brief review on the Kronecker product of two matrix. Denoted by Mm,„(F) 
the set of all to by n matrices over field F, and abbreviate M„.„(F) to M„(F). The Kronecker 
product of A = [aij] G Mm_„(F) and B G Mp.g(F) is denoted by A O B and is defined to be 
the block matrix 



The Kronecker products has many properties^^^l. Here we mention two properties, one of 
which is as follows: 

Theorem 2 (mixed-product) Let A G M„,„(F), B G Mp,g(F), C G M„,fc(F), and D G 
M,,r(IF). T/ien (A ® B)(C (g) D) = AC ® BD. 

Another property is concerned with the eigenvalues of the Kronecker product of two square 
complex matrices. 

Theorem 3 Let A G M„, and B G M,,,,. Denote by ct(A) the all eigenvalues of matrix A. // 
A G cr(A) and x G C" is a corresponding eigenvector of A, and if fJ- (z cr(B) and y G C™ is 
a corresponding eigenvector of B, then A/i G cr(A B) and x (x) y G C"™ is a corresponding 
eigenvector o/A®B. Every eigenvalues o/A(g)B arise as such a product of eigenvalues of A 
and B. 



follows [241: 




From theorem |3] and the fact that the determinant of any matrix is equal to the product of 
its eigenvalues, it follows thafl'^^l. 
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Theorem 4 Let A e M„ and B e M„ be given. Then det(A ® B) = (det A)'"(det B)". 
Thus, A ® B is nonsingular if and only if both A and B are nonsingular. 

In this paper, we need to solve an equation whose coefficient matrix is the Kronecker product. 
The following theorem shows us how to solve the problem'^^' . 

Theorem 5 LetW denote afield. Matrices A e M,„,„(F), B e Mg,p(F), and C e Mm,q(F) are 
given and assume X G M„.p(F) to be unknown. With matrix "K, associate the vector vec(X) G F 
defined by 



vec(X) = [an, • • • a„i, ai2, • • • , a„2, • ■ ■ , flip, • • • , a. 



Then, the following equation: 

(B ® A)vec(X) = vec(C) (3) 

is equivalent to matrix equation: 

AXB^ = C. (4) 
Obviously, equation (j4]) is equivalent to the system of equations 



AY = C 
BX^ = Y^. 



(5) 



3 Algorithms 



3.1 Univariate Interpolation Polynomial 

Let f{x) = X]"=o '^i'^* univariate rational polynomial with degree n. Its approximate 

support points are {{xi, fi)} for i = 0, 1, • • • , rt, which means that 

fi^i) = fi, for i = 0, 1, • • ■ ,n 

where Xi are called interpolation nodes and values fi interpolation datum. We can construct 
polynomial f{x) from its support points by interpolation method. Polynomial interpolation is 
a classical numerical method. It is studied very well for univariate polynomials. In general, 
interpolation problem is essentially to solve the following equation: 



Xq 
Xi 



ai 
\ an J 



= U„+i(xo, ■ • • ,Xn)a : 



/l 



(6) 



where 



U„+i(xo, • ■ • , Xn) 



( 



\ 1 



X3 



xi 



Matrix U„+i (xq, ■ • ■ ,x„ 
detU„+i(a;o, • ■ • ,a;„) = 



a = 



a2 
as 

V a„ y 



(7) 



) is called Vandermonde matrix and its determinant 'Vn+i {xq, ■ ■ ■ ,Xn) — 
Y[o<i<j<ni^j ~ ^i) called Vandermonde determinant. Let D„(z,_;/) 
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denote the determinant of the submatrix of matrix U„+i(a;o, ■ ■ ■ ,Xn) resulting from deletion 
of row i + 1 and column J + 1. D„(i, j) is often called generalized Vandermonde determinant. 
How do we compute it? The follow theorem shows us do it^^^l. 

Theorem 6 Let 5k{xi,X2, ■ ■ ■ ,x„) ^Y,i<ii<i2<---<ik<n^ii^i2 ■'■^ik elementary symmetric 
polynomial with degree k. Then 

From theorem ini it follows that the solution of equation ([6]) 



V„+i(a;o,a;i, • ■ • , a;„) 



(8) 



Theoretically, we can compute the interpolation polynomial after choosing n + 1 distinct 
points , xi , • • • ,Xn and then obtaining their corresponding exact interpolation values /o , /i 7 • • • , fn- 
However, in practice, we often get the approximate function values of f{x), denoted by /o, /i, • • • , /«. 
So an approximate interpolation polynomial f{x) = X]J=o '^3^'^ only produced. The following 
theorem gives an error estimation of coefficients of the approximate interpolation polynomial. 

Theorem 7 Let e = max{|/j — fj\,j — 0, • • ■ ,n}, Xx — min{|a;j — Xi\,i ^ j} and M = 
max{l,max{|a;j|, j = 0, ■ • ■ ,n}}. 

max{|a,-5,|,j-0,--- ,«}<— (n + l)(^^^^2jM", (9) 

where \x\ stands for the greatest integer which is less than or equal to x, and (|^„"2j) number 
of \n/2\- combinations of an n-element set. 

Before giving the proof of theorem [71 we introduce two lemmas: 

Lemma 1 Let 5k{xi,X2,--- ,2;„) = Y,l<ll<l2<■■■<lk<n^^l^^■i ' " elementary symmetric 

polynomial with degree k. Then 

ji n points n+l points 

5k{xo,xi, ■ ■ ■ ,Xj^i,Xj+i, ■■■ ,Xn) = {n + 1- k)Skixo,xi, ■■■ ,Xj,--- (10) 

3=0 

Proof: Every term a; j(, a; • • ■Xi^_^ oi5kixo,Xi, - ■ ■ ,Xj,--- ,x„) does not appear in 4(a:o7 a;i7 • • • , Xi-i,Xi+i,- ■■ ,Xn) 
with i = ii for ? = 0, 1, • • • ,k—l, and appears in the other cases. Therefore, term Xi^Xi-^ ■ ■ ■ Xi^_-^ 
appears {n + 1 — k) times on the left hand side of equation (|10p . So, the lemma is finished. 
Another lemma is concerned with binomial coefficients as follows t^^l; 

Lemma 2 For n a positive integer, the largest of the binomial coefficients 



fdj \IJ \2, 

is (^„"2j)' ^'^^''s L^J stands for the greatest integer which is less than or equal to x. 
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Now we turn to give the proof of theorem [71 
Proof: From equation ([8]), we have 



a, — a 



ELo(-i)'+^D"(i,i)(/i-/.: 



< 



E 



V„+i(xo,Xl, • • ■ ,Xn) 

D„(*,j)ii(/.-/.: 



V„+i(xo,Xi, • • ■ ,Xr. 



|D«(«,j) 



V„+i(xo,xi 



< e 



E- 

i=0 



|V„+i(a;o,a;i, • • ■ ,x„)\ 

n values 



i=0 



noticing equation pO)) yields 



n+l values 



'j\ < ■j^{n+l-n + j)Sn-j{\xQ\,---,\xi^i\,\xi\,\xi+i\, 



\Xn\) 



From lemma m it follows that 



< 



< 



e n! 
AS^''^ ^n-K2j)!Ln/2j! 
e n! 



A' 



n 
n/2\ 



-{n+l){ , 'L, lAf". 



The proof is completed. 

And now, let us discuss how to recover the exact polynomial from its approximate interpo- 
lation polynomial. 

Let n denote the degree of f{x) and N an upper bound of absolute values of denominators 
of its coefficients. Note that f{x) is unknown, so we should estimate an upper bound N and its 
degree n in advance. We may calculate them by the form of the given expression of polynomial 
f{x), and should obtain as less bound as possible. The less bound we obtain, the less the 
amount of computation is for obtaining approximate interpolation polynomial. Once an upper 
bound N and n are gotten, we choose n + l interpolate nodes xo,xi, ■ ■ ■ , a;„ and calculate 



(11) 



Then, compute the values fi w f{xi) for i — 0, 1, • • • ,n with an error less than e. By inter- 
polation method, we compute the approximate interpolation polynomial f{x) with coefficient 
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error less than 1/(2A^^). Finally, use algorithm [T] to obtain the exact polynomial f{x) from its 
approximate polynomial f{x). In summary, the algorithm is as follows: 

Algorithm 2 Input: an expression f{x) which result is a polynomial; 
Output: a polynomial g{x) such that g{x) — f{x). 

Step 1: By the structure of the expression, estimate an upper hound on the degree 
of f{x) and an upper bound on absolute values of denominators of its coefficients, 
denoted by n and N respectively; 

Step 2: Choose distinct points xo,xi,--- ,a;„; Compute = min{|a;j — Xi\,i ^ j} 
and M = max{l, max{|a;i|, i — 0, ■ ■ ■ , n}}. 

Step 3: Compute e in formula hll\) : 

Step 4-' By numerical method, approximately compute the values of f{x) at the points 
Xq,Xi, ■ ■ ■ ,Xn with an error less than e and denote the corresponding values by 
ft ~ f{xi), for i ^ 0,1,- ■■ ,n ; 

Step 5: By interpolate method, obtain approximate interpolation polynomial f{x) ; 

Step 6: Call algorithm[l\to recover the exact coefficients from the coefficients of f{x) 
one by one. Denote the exact polynomial by g{x); 

Step 7: return g{x). 

The correctness of the algorithm above is shown as follows: From step 1 to step 4, we obtain 
the interpolate datum fi with an error less than e. By theorem [71 The errors of the coefficients 
of the approximate polynomial are all less than 1/(27V^), which ensure the exact coefficients of 
the polynomial to obtain from its approximation by algorithm [T] 

3.2 Multivariate Interpolation Polynomial 

For simplicity, we first consider the bivariate interpolate problem, and then generalize the 
results to the case of multivariate interpolation. Let f{x,y) — J2i j '^ij^^V'' be a polynomial 
with rational coefficients, and n, m be the bounds on the degree of f{x, y) in x, y respectively. 
We choose the interpolation nodes {xi, yj) (« = 0, ■ • • , n; j = 0, • • • , m), and obtain the values 
of f{x, y), denoted by £ M (i = 0, • • ■ ,n; j — 0, ■ ■ ■ , m). The set of monomials is ordered as 
follows: 

{xiyj\i = 0, • • • ,n; j = 0, • • • , m} 

{1, y, • • • , y", x, xy, • • • , xy"\ • • • , x", x"y, ■■■ , cc^y™}, 

and the interpolation nodes in the corresponding order is as follows: 

{{x^,yj)\i = 0, • • ■ ,n; j = 0, • • ■ ,rn} = {{xo,yo), {xo,yi),-- ■ , (xo,y™), 
{xi,yo),{xi,yi), - ■ ■ ,{xi,y,n),- ■ ■ ,{xn,yo),--- ,(x„,y„i)}. 

Note that the above ordering is reverse to conventional lexicographic order. Of course it is a 
lexicographic order under assumption of ordering > 1 > 2 > ■ • • and y > x. They do not 
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coincide with our convention, so the above ordering is called reverse lexicographic order in this 
paper. Let 



A = 



00 



fnO 



fon 

fnr, 



the bivariate interpolate problem can be expressed as to solve the following equation: 

Mvec(A'^) = vec(F'^), (12) 
where the coefficients matrix M — 15 x ® Uj^. U^; and Uj, are Vandernionde matrix, i.e., 



/I Xq xl 
1 Xl x\ 



\ I Xn 



J 



/I 2/0 yl ■■■ yJT \ 



1 yi yi 



V 1 y« 



yi 



yZ } 



As said in section 2, in practice we often only get the approximate interpolation values 
/ij ~ fi^iTVj)- Of course, the interpolate polynomial obtained is an approximation. The 
following theorem gives an estimation error: 

Theorem 8 Set fij — f{xi, yj) and denote by fij the approximation of fij. f{x, y) — 'Yliij o.ijX^y'' 
is the interpolate polynomial from interpolate datum fij. Let e = max{|/y — fij\ ■ i ~ 
O,--- ,n;j = 0, •■• ,m}, = min{\xj - Xi\ : i ^ j}, Xy = min{|yj - yi\ : i ^ j} and 
Mx — maxjl, max{|a;j | : j — Q, - ■ ■ , ri}}, My — max{l, max{|yj | : j = 0, • ■ • , Then 



max I Oij — ciij I < 



(m + l)(n + l)M^M" 



A^Ay 



Proof: Let A = (fly) and F = {fij). From equation it holds that 

Mvec((A - A)^) = vec((F - F)^), 
by theorem [SI the above equation is equivalent to the following equation: 



Thus, it is equivalent to 



\]y{A~A)'^\]J = (F-F)^ 

JJyZ = (F F)"^ 
Ux(A - A) = 

Where Z = (zy ). Matrix equation (|13aP is equivalent to 



UyZ., = (F,. -FO^ 



I = 1, • • • m + 1 



(13a) 
(13b) 



(14) 



where Z.^ stands for the i-th column of Z and F;. the i-th row of matrix F. By theorem [7l for 
every i, it holds that 



max \ zji\ < 



Xl 



-(to + 1) 



TO 

Lto/2J 



Ml: 
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Hence, we conclude that 



V 



Let 

argue equation (|13bp in the same way as do above, we deduce that 



m 

"^T'^'" A^A™ \[n/2\J\[m/2\ 

The proof is finished. 

Now, After having a relation between error of approximate coefficients and that of interpolate 
datum, we establish an algorithm to obtain an exact polynomial from its approximation as 
follows: 

Algorithm 3 Input: an expression f{x,y) which result is a polynomial; 
Output: a polynomial g(x, y) such that g{x, y) — f{x, y). 

Step 1: By the structure of the expression, estimate an upper bounds on the degrees 
of f{x,y) in x and y, denoted by n and m respectively. Estimate an upper bound 
on absolute values of denominators of its coefficients, denoted by N ; 

Step 2: Choose interpolate nodes {xi,yj),(i ~ 0,1,--- ,n;j = 0,1,--- ,mj. Compute 
Xx = miiij^i{\xj - Xi\}, Xy = minj^i{\yj - yi\} and = max{l,max"^Q \xi\}, 
My = max{l,max™Q \yi\}. 



Step 3: Compute e: 



\ n \ m 



2(n + l)(m + 1)(l„%j) (L™':;2j)A^."M™iV^ ' 

Step 4-' By numerical method, compute the approximate values of f{x,y) at the points 
{xi,yj)(i — 0,1,--- ,n;j — 0,1,--- ,m) with an error less than e, and denote the 
corresponding values by fij f{xi,yj), (i — 0,1, ■ ■ ■ ,n; j — ■ ■ ■ ,m) ; 

Step 5: By interpolate method, obtain approximate interpolation polynomial f{x,y); 

Step 6: Call algorithm [I] to recover the exact coefficients from the coefficients of 
f{x,y) one by one. Denote the exact polynomial by g{x,y); 

Step 7; return g{x,y). 

As for generalization of the above results to the case of multivariate polynomials r > 2, we 
assert that the situation is completely analogous to the bivariate case. 

Let f{Xi, ■ ■ ■ ,Xr) = '^Oi-^.-.i^Xl^ ■ ■ ■ Xp' be a polynomial in Xi, ■ ■ ■ ,Xr, and let n,; be a 
bound on the degree of f{Xi, ■ ■ ■ ,Xr) in Xi {i — 1, ■ ■ ■ ,r ). Denotes by N an upper bound 
on absolute values of denominators of coefficients of polynomial f{Xi,--- ,Xr). Choose in- 
terpolation nodes {xu-^ , X2i2 ■ ■ ■ , Xri,.) (ife = 0, 1, • ■ • , n^), where k = 1, ■ ■ ■ ,r. Set fi^i2---ir = 
f{xii-^, X2i2, ■ ■ • i^^rv) and denote by /iiia-Jr the approximation of fi^i^-.-i,., by which we com- 
pute an approximate interpolation polynomial f{Xi, X2, ■ ■ ■ , A,.) — ^ Oi-^i^.-.i^X^^ Xl^ ■ ■ ■ X^r . 
Thus we have the following theorem: 
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Theorem 9 Let e ^ niaxjj^i2_...i^ |/iii2---v - f^li2■■■i,^J = min,j^j \xkj - Xki\, and Mk = 
maxj 1 , niax"^Q | Sfc j | } . Then 

max |aiii,2i2,---,r-i, - aiii,2i2,--,rtj < e tttttt • (15) 

il,^2,■■■.^. ^J-^ (Aj^, ) 

Proof: By the Reverse Lexicographic order, the interpolation problem comes down to solving 
the following equation: 

(Uxi «) ■ ■ ■ «) Uxja = F, 

where a is a column vector consisting of elements fliiia - ir reverse lexicographic order and 

F is a column vector consisting of elements fi^i^-.-i^ in the reverse lexicographic order. Thus we 
have the following equation: 

(Uxi ® ■ • • ® UxJ(a - a) = (F - F), (16) 

where a and F are column vectors respectively consisting of elements ai^i^...i^, and fi^i^-^-i,. in 
the reverse lexicographic order. Equation (fTHl) is equivalent to the following equation: 



((Uxi ®Ux2 •■•®Ux._J<5§UxJ(a-a) - (F - F). 

Hence, by recursion and as did in the proof of theorem [51 we can deduce that inequality (|15p 
holds. The proof is finished. 

Based on the above theorem, algorithm [3] is generalized to the case of multivariate interpo- 
lation as follows: 

Algorithm 4 Input: an expression f{Xi, ■ ■ ■ ,Xr) which result is a polynomial; 
Output: a polynomial g{Xi, • • • , Xr) such that g{Xi, • ■ • , Xr) = f{Xi, ■ • ■ , Xr). 

Step 1: By the structure of the expression, estimate bounds on the degrees of f{Xi, • • • , Xr) 
in Xi , denoted by Uifi — 1, ■ ■ ■ ,r) . Estimate an upper bound on absolute values of 
denominators of its coefficients, denoted by N ; 

Step 2: Choose interpolation nodes {xu^, ■ ■ ■ ,Xri^),(ik = 0, 1, • • • ,nk;k = 0, 1, • • • ,r). 
Compute Afc = minj^i{\xkj —Xki\} and Mk — max{l,max"^g \xki\} fork = 1, 



, r. 



Step 3: Compute e: 

nLi(Ar) 



2^^nLiK + i)^r(L„r/2j) 



(17) 



Step 4-' By numerical method, compute the approximate values of f{Xi,--- ,Xr) 
at the points (xu-^,--- ,Xri,.)(ik — 0, 1,-- - ,nk',k = 0, 1,-- - ,r) with an error less 
than e and denote the corresponding values by fi-^-.-i,. ~ /(xii^,--- ,Xri^), fik ~ 
0,1,--- ,nk;k = 0,1,- ■■ ,r) ; 

Step 5: By interpolate method, obtain approximate interpolation polynomial f{Xi, • • • , Xr); 

Step 6: Call algorithm [I] to recover the exact coefficients from the coefficients of 
f{Xi, • • • , Xr) one by one. Denote the exact polynomial by g{Xi, • ■ • , Xr); 

Step 7: return g{Xi, • • ■ , Xr). 
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The correctness of the above algorithm is ensured by theorem [8] and theorem [T] The proof 
can be given similarly to that in algorithm [3l 

4 Experimental Results 

Our algorithms are implemented in Maple. The following examples run in the platform of 
Maplell and PIV3.0G,512M RAW. The first three examples illuminate how to obtain exact 
interpolation polynomials. Example 4 tests our algorithm for more variables and higher degree 
interpolation multivariate polynomials. 

Example 1 Let f{x) be an unknown univariate rational polynomial. Assume that we know 
in advance the degree of f{x): n — 8 and an upper bound of absolute values of denominators of 
its coefficients N — 181. According to theorerr{71 computing exact interpolation polynomial f{x) 
as follows: Choose floating-point interpolation nodes[0.5001, 1.5003, 3.1201, 2.314, 4.02, 5.23, 6, 6.8, 
7.2001]; Calculate A^^ = 0.4001, M=7.2001, Compute e = 0.2202 x 10~". Compute the ap- 
proximate interpolate datum /,; « f{xi) such that \fixi) — fi\ < £, for i = 0, • ■ • , 8, and use 
interpolation method to obtain approximate interpolation polynomial f{x): 

.33333333333333327598X* - .027777777777775996307a;^ -I- .66666666666664376920x^ 
-.072222222222064098893x5 + 1.9583333333326967612a;'' -f 1.7500000000015123646a;^ 
-.24166666666870050587x2 -t- .17500000000137107868a; -f 1.4999999999996690437. 

Then we recover the exact coefficients of interpolation polynomial by algorithm [TJ 

^, , 1 g 1 7 2 g 13 5 47 4 7 3 29 2 7 3 

ffx = -a;** x^ + -a;^ a;^ H x^ + -x"^ a;^ H x + -. 

^ 3 36 3 180 24 4 120 40 2 

Example 2 Let f{x,y) be an unknown bivariate rational polynomial. Assume that we 
know in advance the degrees of f{x,y) in a; and y: n — 3,m = 3 , and an upper bound of 
absolute values of denominators of its coefficients N = 13. According to theorem [51 computing 
exact interpolation polynomial f{x,y) as follows: Choose floating-point interpolation nodes 
X = [0.1,0.5,1.2,1.3], Y = [0.2,0.8,2.1,2.6], and calculate = 0.1, Aj^ = 0.5,14 = 1.3, = 
2.6; Compute e = 0.6659 x 10~^°. Compute the approximate interpolate datum such that 
\f{xi,yj) — fij\ < £, for ? = 0, • ■ • , 3, j = 0, • ■ • , 3, and use interpolation method to obtain 
approximate interpolation polynomial f{x,y) as follows: 

-.2500000301a;3?;2 _ .25000001 lla;^^^ + .08333326696a;2y - .08333336318a;?/2 
-.1666666674y^ - .2499999755a;2 + .2500000031?/^ - .5000000114a; - .5000000029y. 

Then we recover the exact coefficients of interpolation polynomial by algorithm [TJ 

fix.y) — — x y X y -\ xy xy y x H — y x y. 

j\ 4 4 12 12 6 4 4 2 2 

Example 3 Let f(x,y,z) be an unknown tri-variate rational polynomial. Assume that 
we know in advance the degrees of f{x,y,z) in x,y,z: n = 3,to = 2,Z = 3, and an upper 
bound of absolute values of denominators of its coefficients N = 231. According to theorem 
[9l computing exact interpolation polynomial f{x,y,z) as follows: Choose floating-point in- 
terpolation nodes X = [0.5, 1.2, 1.8, 2.7], F = [0.3, 1.5, 2.2], Z = [0.6,1.8,2,2.8], and calculate 
A^ = 0.6, Ay = 0.7, A;, = 0.4, = 2.7, = 2.2, = 2.8; Compute e = 0.735 x 10^1°. 
Compute the approximate interpolate datum fij^. such that \f{xi,yj,Zk) — fijk\ < for 



12 



YONG FENG ■ JINGZHONG ZHANG ■ XIAOLIN QIN ■ XUN YUAN 



i = 0, • ■ • , 3, j — 0, • • • , 2, fc = 0, • • • ,3, and use interpolation method to obtain approximate 
interpolation polynomial f{x, y, z) as follows; 

-.08333331999a;y + .1666668454x^2/22 - .00A62M85562xy^z^ - .02777770742a;yz3 
-.05555579922xy2z - .3333337487a;yz2 + .08333353340x^2/ - .1666659870x2/z 
+.013888668362/Z+ .083333275412^ + .5000001136X+ .1666667565y. 

Then we recover the exact coefficients of interpolation polynomial by algorithm [TJ 

rr \ 132,122 1 22 1 3 I2 

f x, y, z) — X V H — X V z xii z xiiz xy z 

j\ ,y, ; 12 ^ 6 ^ 216 36 18 ^ 

1 2 1 2 1 1 1 2 1 1 

— xyz H X y xyz H yz H z H — x H — u. 

3 i/ -r ^ 6 72^ 12 2 6^ 

Example 4 This example tests algorithm U by computing determinants, and the results 
are shown in table 1, where '?' represents that the computation is not finished. Table 1 shows 
that our algorithm is more efficient than exact interpolation algorithm when scale of problem 
is larger. 



Table 1: Comparison between our algorithm and exact interpolation 



matrix 
sizes 


variates 


terms 


error control 


time(sec) 


exact in- 
terp(sec) 


3x3 


3 


210 


0.1679e- 26 


1.766 


1.531 


3x3 


4 


343 


0.1829e- 30 


5.578 


3.837 


4x4 


4 


256 


0.5441e- 19 


2.078 


1.734 


5x5 


3 


816 


0.4471e - 37 


5.837 


4.241 


6x6 


3 


2774 


0.3781e-45 


147.452 


67.731 


7x7 


4 


5657 


0.4731e - 50 


258.129 


347.372 


8x8 


4 


8847 


0.4251e- 55 


657.221 


899.372 


9x9 


3 


15367 


0.5882e - 63 


2074.771 


2744.305 


10 X 10 


4 


23371 


0.7321e- 61 


5721.491 


8737.707 


11 X 11 


4 


40721 


0.6781e- 73 


9287.527 


? 



5 Conclusions 

The paper presents a new method for obtaining exact interpolation multivariate polynomial 
with rational coefficients from its approximate interpolation. The exact results can be obtained 
by our algorithm as long as we estimate the upper bounds on the degrees in every variable 
and upper bound on absolute values of denominators of its coefhcients. As we know, many 
complicate polynomials computation can be completed by exact interpolation method, so they 
can also be completed by our algorithm, i.e., by approximately numerical computation. 

The paper proposes a way to get an exact polynomial by approximate numerical compu- 
tation. Like exact interpolation methods, this algorithm is a parallel method. The accuracy 
control e in formula (|17p is an exponential function in degrees of polynomial, and is a polyno- 
mial function in the upper bound 1/A^. Experimental results show that our method is more 
efficient than the exact interpolation methods when scale of a problem gets larger. 
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